Parameter estimation from gravitational waves generated by non-spinning binary 
black holes with laser interferometers: beyond the Fisher information. 
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In this paper we apply to gravitational waves from non-spinning binary systems a recently intro- 
duced frequentist methodology to calculate analytically the error for a maximum likelihood estimate 
(MLE) of physical parameters. While existing literature focuses on using the Cramer Rao Lower 
bound (CRLB) and Monte Carlo simulations, we use a power expansion of the bias and covariance 
in inverse powers of the signal to noise ratio. The use of higher order derivatives of the likelihood 
function in the expansions makes the prediction also sensitive to the secondary lobes of the MLE 
probability distribution. We discuss conditions for validity of the CRLB and predict new features 
in regions of the parameter space currently not explored. For example, we see how the bias can 
become the most important contributor to the parameters' errors for high mass systems (200A/q 
and above). 



I. INTRODUCTION 



Coalescing black holes binaries (BBH) are among 
the most promising sources of gravitational waves 
(GW) transients Observations strongly indi- 

cate the existence of stellar-mass BHs (3-10 Mq ) 
H and super massive black holes (10 4 — 1O 1O M0) 
also suggesting the possibility of interme- 
diate BHs (0|)- Observations and models also 
point to the formation of binary BH systems (BBH) 
(1,0,0,0,0,0,0,0,0)- The chances of 
observation of BBH GWs in the next few years are 
promising given the advanced generation of laser inter- 
ferometers currently under construction (0,0)) and 
recent advances in modeling the BBH waveforms. In 
fact, even if the inspiral and ringdown phases of the 
life of a binary system are well understood (the in- 
spiral GW can be computed using Post-Newtonian 
approximations ([HI), while the GW for the ring- 
down phase can be obtained with black hole perturba- 
tion theories) only recently numerical breakthroughs 
in numerical relativity made compute GWs from the 
mergerphase possible (0, 0, 0, 0, 0, 0, 
0, 0)- The GW community is actively preparing 
for detection and parameter estimation opportunities 
in GWs from BBHs 0] and plans are shaping up 
to explore all the interesting regions of the parameter 
space. P. Ajith and collaborators ([38[) have proposed 
a template bank for the waveform coming from a co- 
alescing binary system, made of non-spinning masses, 
that takes into account the inspiral merger and ring- 
down stages of the binary's life (IMR). IMR wave- 
forms are obtained by tuning numerical parameters 
of a phcnomenological wave (see [38j ) to a set of nu- 
merical calculations performed in the full GR. Other 
attempts to produce analytical approximations have 
been developed in more recent times, like the effec- 
tive one body model (EOB) proposed in 0], that 



was used in the fifth LIGO run S5. The most accu- 
rate EOB model in the literature is described in (0, 
0])- The discrepancies between different analytical 
approximations to numerical relativity waveforms are 
described in 0] and show for example that the ampli- 
tude of the IMR waveforms can be up to 20% different 
from the EOB ones. This means that a binary at a 
given distance could produce a signal to noise ratio 
(SNR) at receiver up to 20% different from the one 
expected using IMR waveforms. 

This paper discusses the estimation accuracy of 
BBH physical parameters which can be obtained with 
advanced configuration laser interferometers and IMR 
waveforms. Existing literature 0] predicts the er- 
rors with approaches that have some intrinsic limi- 
tations: (a) with the square root of the inverse of 
the Fisher information matrix elements (commonly 
named Cramer Rao Lower Bound (CRLB), 0]) that 
only takes into account the curvature of the likelihood 
function around the true value of the parameters. The 
CRLB is known to underestimate the error in low sig- 
nal to noise ratio (SNR) (0, 0, 0) and is the 
lowest possible uncertainty for unbiased estimators, 
or (b) perform simplified Monte Carlo simulations. 
In the MC simulations r] is enforced to values 
<= 0.25 and the simulations do not explore the sec- 
ondary peaks of the likelihood function. 

A recent paper, f22j . from the authors, studies the 
problem of MLE errors from GWs from the inspiral 
phase of binary systems with asymptotic expansions 
for the covariance and the bias in terms of power series 
in the inverse of the SNR. The first order of the co- 
variance series is the inverse of the Fisher information 
matrix. The second order is a more complicated ex- 
pression that depends on the secondary maximum of 
the probability distribution because it contains higher 
order derivatives of the likelihood function (up to the 
fourth). The first two orders of the covariance and 
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bias expansion are a better tool than the CRLB to es- 
timate the errors and allow the determination of nec- 
essary conditions for the validity of the CRLB (for 
example by requiring the second order to be much 
smaller than the first). These conditions allow predic- 
tions to be made on the interferometers' capability to 
estimate parameters in different regions of the param- 
eter space. We show that the variation of the errors in 
the parameter space is more complicated than as pre- 
dicted solely by the CRLB and that the existence of 
minima for particular values of the masses can hap- 
pen even without including the merger phase. We 
also predict that the bias can become the most im- 
portant contributor to the parameters' errors for high 
mass systems (200M© and above), due to the nonlin- 
ear dependence of the signal on the parameters (these 
regions of the parameter space are not yet explored) . 

In section [TT| we define the data model, the MLE, 
the statistical errors and how to compute the asymp- 
totic expansions. In section 11111 we give the analytical 
expression for the IMR GW for BBHs, and the noise 
spectra for the advanced configurations of LIGO and 
Virgo. In section [IV] the results are described. 



II. STATISTICAL MODEL 

We model the output x(t) of a GW detector as the 
sum of the GW signal h(t, M ), that depend on the 
vector of unknown parameters 0^ , and a Gaussian sta- 
tionary noise w(t) with zero mean, E[w(t)] = 0, with 
E[A] denoting the mean of A over the ensemble 



of two functions can be written like a scalar product 



x(t) = h{t,6») + w{t) 



(2.1) 



If one expects the wave to have a known analytical 
form it is possible to perform an estimation of the 
parameters by filtering detector data with a bank of 
waveform templates. In Gaussian noise, this is a max- 
imum likelihood estimation. 

Given the Fourier transform for a function h(t) as 
h(f) = J dte~ 2lTl f t h(t) 1 the expectation of the product 
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(«(/) I <f)) = 2 
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df 



u(fHfr+u(frv(f) 

s h (f) 



(2.2) 

where the range of integration depends on the antenna 
properties and on the theoretical model for the binary 
system, and where we introduced the one sided noise 
spectral density, Sh(f) defined by E[w(f)w(f')] = 
\S h (f)5{f - /'). The SNR corresponding to the op- 
timal filter is defined as 



W),h(f)) 



(2.3) 



Once the values of the parameters are estimated using 
matched filters, the accuracy can be evaluated with 
the square root of the mean squared error (MSE) for 
the j-th parameter: 



MSE& = E 



E 



■do - 



- ~ 2 _L h 2 

= a ,v + b ,v 



(2.4) 



For large SNRs one can use the CRLB to obtain 
a lower bound for the error of the j-th parame- 
ter MSEtfj > • •, where i is the Fisher In- 
formation matrix, whose (j k) element can be writ- 
ten as a scalar product of signal's first derivatives 
i jk ee (h(f)j,h(f) k ) where h(f), = ^f. The 
above notation allows for more concise formulae: 



(hal 



(ab- 



■ p) , where the f de- 



pendence is not explicitly shown. The scalar product 
is the one defined in cq. 12.21 In [22j expansions for 
both the covariance and the bias like power series on 
1/p are given as: 
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= ^ < [l]+4 i [2] + ---(2.5) 



Ml]+M2] + --- (2.6) 



where the first order covariance is the CRLB, and the 
second order can be written as: 



o% [2] = i 3m i 3n i pq (v nmpq + 3(nq \pm) + 2v nmv<q +v mvq>n ) + 



(2.7) 
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where 



E [ftoi— o, hb ± — b p ■ ■ ■ h Zl ... Zq ] 



are expectations of product of the signal derivatives in 
the time domain. For example: u a &.c.d = E[h a b h c hd\. 
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Using some algebra, the following explicit expressions v ab ,cd,e = —iabVcd , e — icd^ab , e (2-17) 
can be given, in the frequency space: Vabcde = _( abcd | e ) _ ( a6ce | d ) _ ( aMe | c ) _ ( acde | 5) 

— (bcde I a) — (abc \ de) — (abd \ ce) — (acd \ be) 

— (bed I ae) — (abe \ cd) — (ace \ bd) — (bee \ ad) 
-(ade I be) - (bee | ac) - (cde | ba) (2.18) 



V a ,b = - Vab = lab = | &) (2.8) 

v a b c — (ab \ c) (2.9) Where ijk is the Fisher information matrix. 

v a bc,d = (abc | d) (2-10) 

Vabc = -(ab | c) - (ac | b) - (be \ a) (2.11) 



Vabcd 



The reason for the presence of 3(nq \ pm) in (|2.7[) is 
that this expression is a simplified version of the eq. 
2.4 we gave in (22|, in which some of the v have been 



v a b , cd = (a& | cd) + (a\b)(c\ d) (2.12) replaced with their values in terms of scalar product. 



-(ab | cd) - (ac \ bd) - (ad \ be) - This allows for slightly shorter expression. The same 

kind of simplifications can be performed on the first 
and second orders of the bias, that have a final form 
shorter than that presented in [22j: 



— (abc | d) — (abd \ c) — (acd \ b) — (bed \ a) 
v a b, c ,d = -(a\b) (c\d) = -i a bicd ( 2 -13) 
Vabc , de = (abc | de) - i de v abc (2- 14 ) 



v abcd!e = {abcd\e) (2.15) 6*r[l] = \i ra i bc (v abc + 2v c , ab ) (2.19) 



V a bc,d.e = ideV a bc (2.16) 
I 



2 



b^[2] 



[vabcde + 4(ac | bde) + 8(de | abc) + Av abc ^ d ] 



{2v a f ed Vgb,c + 2VbedfV a c,g + 4v a bedVgf,e) + (v a f ed V gch + 



+ 2v ab edVgcf + 2v dbeg v acf ) + (2v aed (gb | fc) + Av acf (dg I eb) + Av bed (ac | gf) 
+ 2v fcb (ag | ed)) + {4v a f e , g Vdb,c + &v a f etC v dbtg + iv db e,gV a f,c) + (2v abe , g v cd f 
+ 4v dbe ,gVacf + ^VabejVcdg + 2v dge . b v acf ) + (A(ag | fc) v edjb + 4(ed | fc) v ag . b 



A(ag | ed) v fc , b ) 

■ma-bc-de^fg-ti 

8 



)adf(VebcV g ti + 2v etc V gbl + 4:Vgb e V t ci + %V gbt V eci + 2v ebc V gt , 



+ 4w etc w gh , J + 2v gti v eb , c + 4:V gtc v ebii + 8v gbt v C£ti + 8v gbt v ci , e + 8v gbe v ct ^ + 8v cte v gb ., 
+ ^v cU v gb ^ e + 4Vgt,iV eb , c + 4v eb:i v gt , c + 8v gtjb v ic , e + 8v gt:e v lC:b + Av bet v g ^) 
+ v dcl (8v bgt v aeJ + 4v bgf v at , tt + 8v ae ^v bgJ + 8v ae jv b g tt + 8v aftb v ge ^)} 



(2.20) 



III. THE IMR WAVEFORM 

IMR waveforms were obtained by tuning numeri- 
cal parameters of a phenomenological wave to a set 
of numerical calculations performed in the full GR, 
|3^ | . Afterwards, the faithfulness of the IMR waves 
has been improved (jHI, [4~ 
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and they have 
been used for the purpose of parameter estimation. In 
this work we follow the model presented in . The 



(while the EOB are calculated in the time domain) , as 
a piecewise function with a part for the inspiral, one 
for the merger, and another for the ringdown phase. 
Explicitcly, we write the GW in the Fourier space as: 



h(f)=A eff (f)e i *°ffW 



(3.1) 



where the phase and the amplitude arc expressed 



IMR waveform is written directly in the Fourier space as: 
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(// fmerg) ^ if 

A eff (f) = A f-Ul \ (/ / f merg y 2 / 3 if /, 

< Cd£(f, fringe <j) if fring< f < f cut 



f ^ fmerg 
rnerg — f ^ fring 



* eff (f) = 2nft a 



k in {0,2,3,4,6} 



(xkV 2 + VkV + Zk) (ttM/) 3 



(3.2) 
(3.3) 



Wo have defined: 



to 



It® ( fring 



■rnerg 



£{f, fring, f) — 



f 



2tt (/-/„„ ff ) 2 + a 2 /4 



The phenomenological parameters f merg , fring, feat 
depends on the total and symmetrized mass, via the 
following expressions: 



fi = 



aiV 2 + biV + Cj 



I 

vLigo one sided noise spectral density is written as: 

1 -x 2 +x 4 /2 



S h (f) = s a 



.t- 4 - 14 - 5x- 2 + 111- 



1 + x 2 /2 

( 3 - 4 ) S h (f) =oo, f<f low 

where the lower frequency cutoff value is fi ow = 10Hz, 
fa = 215Hz, and 5 = lQ-^Hz" 1 . While for 



/ ^ flow 

(3.6) 



fa 



the AdvVirgo: 
(3.5) S h (f) = So 



2.67 1(T 7 x- 5 - 6 + 0.68 e^ 73 V' 34 



+ 59e (ln;r)2 [- 3 - 2 - imlnx -°- 13 ( lnx ) 2 } x - i - 1 + 



+ 0.68e- a73 ( lna; ) 2 x 5 - 34 



f>fh 



S h (/) = oo, f<fi ow (3.7) 

where the lower frequency cutoff value is chosen to 
be f low = 10Hz, x = /o = 720Hz, and S a = 



10- 47 Hz _1 . Fig. [T] shows the value of y^Sjj) for 
both detectors. 



with i = [merg, ring, a, cut] (for the values of the nu- 
merical coefficients a, b, c, x, y, z sec [39]). They repre- 
sent the frequency at which the system passes from 
its inspiral phase to the merger {fmerg), from the 
merger to the ringdown (fring), and the frequency 
for which the signal ceases to be described by this 
model (fcut, this is also the upper limit of the integrals 
(12.21) ). This signal depends on five physical parame- 
ters (A,t a ,4> a ,M,ri): (I) A is the amplitude of the 

wave. It can be expressed as A = ^2/3 yfli where 
d is the effective distance of the binary. (II) t a is the 
arrival time of the GW at the detector. (Ill) <fi a is the 
arrival phase, i.e. the phase of the signal at the time 
t a . (IV) M is the total mass of the binary. (V) r\ is the 
symmetrized mass ratio: r\ = m! m 2 /M 2 . Sometime 
the chirp mass is used in the literature, instead of the 
total mass. They are related by M c hi rp = tjs M. 
If one considers the merger and ringdown phase too, 
the amplitude is not uncoupled from the other pa- 
rameters. One can not work in the simplified four- 
dimensional parameter space obtained by treating the 
amplitude as a known constant ([22|, [43|,[44|), and 
the full five dimensional space must be considered. 

We perform the calculations using cither the de- 
sign advanced Ligo (AdvLIGO) or advanced Virgo 
(Adv Virgo) noises, as they are given in [4(j. The Ad- 
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Figure 1. Power spectral densities for AdvLigo (continuous 
line) and AdvVirgo (dashed line) 

IV. RESULTS 

To be able to compare our results with those of 
HI, we consider 77 = 0.16,0.2222,0.25 and M = 
20, 100, 200M Q . Tables 1 and 2 show the values of 
the first two orders of the covariance and bias for 
SNR=10. For each value of n and M & the error esti- 
mation for t a , 4> a , total mass and symmetrized mass 
ratio rj are presented. 
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Figure 2. The errors in an Advanced Ligo detector (table 
above) and Advanced Virgo (below). a[l] and a[2] are the 
first (the usual CRLB) and second order in the variance 
expansion; while 6[1] and b[2] are the first and second or- 
der of the bias (see 12.51 and I2.6|l . The time errors are in 
milliseconds, the phase errors are in radians, while the er- 
rors in the mass parameters are in percent. The SNR is 
equal to 10 
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The CRLB consistently underestimate the error at 
this SNR. The errors for different values of the SNR 
p can be obtained by multiplying the first orders for 
10/p and the second orders by 100/p 2 (and multiply- 
ing for 10 and 100 gives Sf and A necessary 
condition for the validity of the CRLB can also be 
obtained as: 



Si 



Sf * 



< 1 



(4.1) 



where the < sign can be replaced by << depending 
on the accuracy needs. Using the data in the tables 
above, we can plot the errors against the SNR, for 



fixed values of the total mass M and the mass ratio 
77. For example Fig. [3] shows the first order variance 
for the arrival time estimation, and the total variance 
(first plus second order) for an equal masses system of 
total mass M = 200 A/©. It is clearly visible for which 
SNR the CRLB ceases to be faithful. 

Another useful application is the calculation of the 
CRLB and second order variance for systems having 
mass and mass-ratio within a chosen range, building 
then a grid of results that clearly show for what sys- 
tems the errors are smaller. We have done that for 
M = 4Af ..2OOAf Q and 77 = 0.10..0.25. In Fig. g]we 
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Figure 3. The error in the estimation of the arrival time for 
a system with M = 2OOM0 and r\ = 0.25. The dot-diamonds 
line is the first order variance (CRLB) while the dashed line is 
the total variance. The bias is not shown, being negligible. 
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Figure 4. The percent error on the total mass estimation (first 
plus second order) as a function of the total mass and sym- 
metrized mass ratio. The systems have a fixed SNR, p = 10, 
and the AdvLigo noise is used. This plot does not depend on 
the values of t a and <j> a 



show a contour plot of the total error in the estimation 
of the system's mass (in percent), for systems having 
a fixed SNR of 10 and using the AdvLigo noise. Fig. 
H] shows a not monotonic trend for the error: for ex- 
emple, a system having rj = .25 will have an error of 
about 1% if its mass is 4M Q . As the mass increases, 
keeping -q constant, (this is equivalent to scanning Fig. 
0] from the lower right corner to the upper right one) , 
there will be a local maximum of the error, due to the 
"island" on the right side of the plot, for M ~ 55M , 
where the error reaches about the 5%, then the error 
goes down for a while, and begins to grow again, from 
about M = 100M o . On the other hand, if 77 = 0.16, as 
the system's mass increases it won't find any "island" 
of great error, and no local maxima will be present. 
We will recover this behavior later (Fig. [5]). 

We can use this kind of plots in another way. We 
calculate the ratio between the total variance (first 
plus second order) and the CRLB, and we let the 
SNR vary until the ratio goes under a chosen thresh- 
old for each point in the grid. Let us for exemplc plot 



'/ 



Figure 5. The ratio between the total mass error (first plus 
second order) and the CRLB, as a function of the total mass 
and simmctrized mass ratio. The systems have a fixed SNR, 
p = 10, and the AdvLigo noise is used. 



^ a[1 l\i\M [2?M - for p = 10 > this is shown in Fi §- G3 

For highly symmetric massive systems the CRLB is 
nearly one half of the corrected error. If we want the 
ratio for the estimation of the total mass to be smaller 
than 1.05 for each value of M and rj in the range con- 
sidered, the SNR must be p > 56.4. We can perform 
the same kind of calculations for the other parame- 
ters. The biggest SNR we calculate in this way is the 
one required for the CRLB of the time parameter to 
attain a 5% precision, p = 61.6. Then we can say that 
for every binary system having mass and symmetric 
mass ratio in the range given above, a SNR of 61.6 
assures that when using the CRLB for the errors' es- 
timation our result will not differ more than 5% from 
the corrected errors. 

It must be stressed that plots like those in Fig. [4] 
and Fig. [5] does not depend on the actual value of the 
arrival time and phase. The reason is the particular 
form of the signal, cq. (|3.1[) . and the fact that t a 
and 4> a are only contained linearly in the phase of the 
signal. Let us for exemple consider the element 
of the CRLB, it depends on the real part of hih*. 
Developing the derivations we have: 

$t[hih*} = di [(Ai + A (tyi)) e** (A, + A (-#,•)) e-**] 

(4.2) 



A 2 ^j 



The arrival time and arrival phase will not be con- 
tained in terms like A or Ai, as the amplitude docs 
not depend on them. Neither they will be in terms 
like ip a , because they were contained linearly in ip. It 
is easy to see that the same kind of proof holds while 
calculating the v..., eqs. (|2.8[) to (|2.18|) . and the opti- 
mal SNR, eq. 

As a general trend above SNR=20 the results are 
consistent with the results derived in [4£j with both 
Monte Carlo simulations and the CRLB. Lower SNRs 
result in higher uncertainties, and occasionally the 
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Figure 6. The relative contribution of the different phases to 
the total squared SNR for an equal mass system of total mass 
M 



Figure 7. The relative contribution of the different phases to 
the total squared SNR (low mass detail) 



Monte Carlo simulations are below the CRLB. This 
behavior has already been observed in (37j . for inspi- 
ral signals, where the author shows that the inconsis- 
tencies between Monte Carlo and CRLB are a con- 
sequence of the restriction r\ < 0.25 in the templates 
bank. By incorporating templates with 77 > 0.25, 37| 
obtains a good agreement between Monte Carlo and 
CRLB for 1.4 — 1.4 and 5 — 5 solar masses systems. 
However, discrepancies are still present in 10 — 10 
solar masses systems. The author acknowledges that 
there is not a satisfactory explanation for this incon- 
sistency More efforts must be made to fully under- 
stand how the boundary 77 = 0.25 affects Monte Carlo 
simulations and CRLB. 

A comparison with our results that involves only 
the inspiral phase [22| indicates that by adding the 
merger and ringdown phases the errors decrease, but 
the necessary SNR (defined below) for the covariance 
to attain the CRLB can be up to a factor of two. For 
example looking at Fig. 8 of [22j 3 of the 4 rows are 
comparable (not the third one because in the IMR 
we use the total mass while in [22| the chirp mass). 
The SNR necessary to attain the CRLB in the IMR 
signals at Mq = 20 is between 9 and 11 while for inspi- 
ral phase waveforms is between 4 and 7. The absolute 
values of the errors are however larger for inspiral sig- 
nals. In this case the contribution to the SNR of the 
inspiral phase is 99 percent (see for example Fig. [6] or 
its low Mq blow up in Fig. [7]). 

The bias does not play an important role, except 
for small SNRs (< 10) or for high mass systems, for 
which the error in the total mass and arrival phase 
is seriously affected by the bias. Unfortunately, while 
the use of the IMR allows for smaller values for the 
error estimation of the arrival time, total mass, and 
77, the same cannot be said of the arrival phase, for 
which the inclusion of the merger and ringdown phase 
seems to degrade the estimation, so that the error of 
the arrival phase estimation is in general higher than 
2-7T, indicating that the arrival phase is unpredictable. 
In [44| the error for the arrival phase was estimated 



using the inspiral 3.5 PN wave, obtaining a value of 
Acj) = 1.16 rad for a system of M = 2OM at an SNR 
p = 10 using the AdvLigo noise. For the same system 
using the IMR wave we obtain A<fi = 14.8 rad (9.98 
rad considering the CRLB only). An estimation of 
the arrival phase error was not given in [4p| , and so a 
direct comparison is not possible. Since they are con- 
sistently above 27r, the errors on <fi are not included in 
the plots. Let us stress that the large errors on the 
arrival phase does not imply that we cannot believe 
the estimations we have for the other parameters. The 
reason is that a correlation coefficient different from 
zero, say close to +1, tells us that an overestimation of 
the parameter x comes together with an overestima- 
tion of the parameter y; but it doesn't tell anything 
abut the relative magnitude of the errors. It is inter- 
esting to plot the errors values against the total mass 
of the system, for a fixed value of the SNR. Fig. [9] 
shows these plots for a mass range from 4M© up to 
5OOM0, and an SNR p = 10, using the AdvLigo noise. 
Fig. [10] does the same with the Adv Virgo noise. The 
inclusion of the second order variance and bias has 
visible consequences on the errors for large mass sys- 
tems, for which the corrected error can be much larger 
than the CRLB. The plots show an oscillatory char- 
acter of the bias, due mainly to the behavior of b [2]. 
For example we plot in Fig. [TTJthc bias orders on the 
total mass estimation against the system's mass, for 
77 = 0.2222. 

The second order covariance and the bias also seem 
to reinforce the minimum of the errors for M 100M Q . 
It is important to notice that a local minimum would 
be present also with the inspiral phase only. For ex- 
ample in Fig. [8] we show the error in the timing calcu- 
lated for a 1.4-1.4 Mq binary system where only the 
insp iral phase is used (3.5 Post-Newtonian waveform, 
[44|). and a minimum is visible once the second or- 
der is included. To verify our new predictions in the 
errors for very large masses, numerical simulations or 
direct applications of parameter estimation pipelines 
are needed. Fig. 9 of [4(| shows how BHs merger 
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Figure 8. Minima in the error with respect to the mass can 
appear because of the contribution of the second order. In the 
plot the timing error is evaluated only for the inspiral signal of 
■q = 0.25 at a fixed SNR, p = 10. 



can have very high SNRs, even of 100, at 100 mcga- 
parsccs for advanced LIGO. At these SNRs our cor- 
rections would not be important. However, in order 
to have useful detection rates we would need to rely 
on sources up to a gigaparsec. In fact, in [45| a re- 
alistic rate for BH merger is given as 0.4 per million 
years in an equivalent Milky Way galaxy. Equation 
5 in (45j also gives an approximation of the number 
of equivalent Milky way galaxies. Combining these 
two observations suggests that a realistic rate for BH 
mergers within a gigaparsec is about 1.6 per year. 
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Figure 9. (Color Online) The errors plotted against the 
total mass, using the AdvLigo noise, for a fixed value of 
p = 10, and with r? = 0.25 (top), 77 = 0.2222 (middle) 
and rj = 0.16 (bottom). The "tot" dot-dashed line repre- 
sents the MSE (variance plus bias); the "var" dashed line 
the variance (first plus second order); the dot-squares line 
is the fist order variance (CRLB); the dot-diamonds line, 
finally, is the absolute value of the bias (first plus second 
order). The "tot" and "var" lines are nearly superimposed, 
except for very high masses 




V. CONCLUSIONS AND FUTURE WORK 

In this work we apply a recently derived method- 
ology to the errors in estimating physical parameters 
from IMR signals. The asymptotic expansion of the 
bias and the covariance are critical to have realistic es- 
timates of the error for SNR below 20 where the first 
detections of present and future laser interferometers 
might live. The behaviour of the errors, in terms of 
minima and maxima, in different regions of the param- 
eter space appear to be more elaborate than predicted 
by the CRLB. For example the bias can become dom- 
inant for MLEs on sy terns with large masses. This 
paper will aid the preparatory work that the scientific 
community is undertaking to prepare for the scientific 
runs of the advanced version of the earth based laser 
interferometers. 
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Figure 11. (Color Online) The bias, first order (red dashed), 
second order (cyan dotted) and total (blue line) using the Ad- 
vLigo noise, for a fixed value of p = 10, and with r\ = 0.2222 
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